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Abstract 


In the current study, a method based on Fourier analysis is 
developed to analyze the force and moment data obtained in large 
amplitude forced oscillation tests at high angles of attack. The aerodynam- 
ic models for normal force, lift, drag and pitching moment coefficients are 
built up from a set of aerodynamic responses to harmonic motions at 
different frequencies. Based on the aerodynamic models of harmonic data, 
the indicial responses are formed. The final expressions for the models 
involve time integrals of the indicial type advocated by Tobak and Schiff. 
Results from linear two- and three- dimensional unsteady aerodynamic 
theories as well as test data for a 70-deg delta wing are used to verify the 
models. It is shown that the present modeling method is accurate in 
producing the aerodynamic responses to harmonic motions and the ramp 
type motions. The model also produces correct trend for a 70-deg delta 
wing in harmonic motion with different mean angles-of-attack. However, 
the current model cannot be used to extrapolate data to higher angles-of- 
attack than that of the harmonic motions which form the aerodynamic 
model. For linear ramp motions, a special method is used to calculat the 
corresponding frequency and phase angle at a given time. The calculated 
results from modeling show higher lift peak for linear ramp motion than for 


1 



harmonic ramp motion. The current model also shows resonably good 
results for the lift responses at different mean angles of attack. 
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Nomenclature 


A. 

Bj 



j 

k 

M 


coefficient of cosine Fourier series 
coefficient of sine Fourier series 

the average value of the constant terms in the harmonic 
oscillation responses 
reference values 

2- D lift coefficient 

3- D drag coefficient 
3-D lift coefficient 

variation of lift coefficient with respect to angle of attack 
variation of lift coefficient with respect to pitch rate 
3-D pitching moment coefficient 
3-D normal force coefficient 
a constant 

constants associated with the virtual mass effect 
constants in amplitude function to be determined 
imaginary part of a complex number 
index 

reduced frequency (=0l>C/v oo ) 

Mach number 
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n index for reduced frequency, also index for the coefficient 

in Padti approximants 
N the number of frequency 

PDj Pad6 approximants 

Py coefficients for Pad6 approximants 

t time 

t’ nondimensional time ^tv^C) 

UQVLM unsteady quasi-vortex lattice method program 

v M free stream velocity 

Greek: 

a angle of attack (=<*0 coskt’) 

defined as c^+a 

oCq amplitude of angle of attack 

a*, mean angle of attack 

tii time rate of change in angle of attack 

a time rate of change in tic 

8 reference length 

t dummy time integration variable 

£ running variable in time 

0 defined as 0=kt’ 

<|> phase angle 
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Chapter 1 
Introduction 


Due to the requirement of increased performance and 
maneuverability, the flight envelope of a modem fighter is frequently 
extended to the high angle-of-attack regime. Vehicles maneuvering in this 
regime are subjected to nonlinear aerodynamic loads. The nonlinearities 
are due mainly to three-dimensional separated flow and concentrated vortex 
flow that occur at large angles of attack. Accurate prediction of these 
nonlinear airloads is of great importance in the analysis of a vehicle’s flight 
motion and in the design of its flight control system. As Tobak and Schiff 
mentioned in ref. 1, the main difficulty in determining the relationship 
between the instantaneous aerodynamic load on a maneuvering vehicle and 
the motion variables is that this relationship is determined not only by the 
instantaneous values of motion variables but also by all of the prior states 
of the motion up to the current state. Due to advanced computer 
techniques, one straightforward way is to solve the flow-field problem and 
the dynamic equation together. For example, a CFD method can be used 
to solve the Navier-Stokes equations governing the separated flow field. 
Then the calculated forces and moments are used in the dynamic equations 
governing the vehicle’s motion to calculate motion variables. The motion 
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variables will change the vehicle’s attitude, and thus the forces and 
moments. Results of repeatedly calculating these coupled equations would 
be the complete time histories of the aerodynamic response and of the 
vehicle’s motion. Although solving these coupled equations is the exact way 
to account for the time-history effects in predicting the aerodynamic 
response to arbitrary maneuvers, this is obviously a very costly approach. 
In particular, at high angles of attack, the aerodynamic loads depend 
nonlinearly on the motion variables. Under such conditions, even if the 
vehicles start from closely similar initial conditions, they can experience 
widely varying motion histories. Thus, a satisfactory evaluation of the 
performance envelope of the aircraft may require a large number of coupled 
computations, one for each change in initial conditions. Further, since the 
motion and the aerodynamic response are linked together in this approach, 
there can be no reutilization of the previously obtained aerodynamic 
reactions. 

To avoid the disadvantage of solving the coupled flow- field equations 
and aircraft’s motion equations, an alternate approach is to use a 
mathematical modeling to describe the steady and unsteady aerodynamics 
for the aircraft’s equations of motion. Ideally, with a mathematical model, 
an evaluation of the aerodynamic terms specified by the model would be 
required only once. The specified model can be reutilized to solve the 
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aircraft’s equations of motion over a range of motion variables and flight 
conditions. 

In the classical linear potential flow theory (refs. 2 and 3), 
researchers in the field of aeroelasticity used the Fourier transform to 
relate the aerodynamic response of step change in angle of attack of a wing 
to that of harmonic oscillatory motions. The transient aerodynamic 
reaction to a step change is termed the indicial function and has been 
calculated for several classes of isolated wings (refs. 2-5). By a suitable 
superposition (ref. 6) of these results, the aerodynamic forces and moments 
induced in any maneuvers can be studied (refs. 2 and 3). Tobak has 
applied the indicial function concept to analyze the motions of wings and 
wing-tail combinations (ref. 7). Later, based on a consideration of function, 
Tobak and his colleagues (refs. 1 and 8) have extended the concept of 
indicial function into the nonlinear aerodynamic regimes. The simplest 
nonlinear aerodynamic model proposed in ref.l has been applied by several 
authors (refs. 9-13) to perform the analysis. However, that simplest model 
is accurate only to the first order of frequency. It needs to be improved for 
a more general response. 

Aerodynamic forces and moments acting on a rapidly maneuvering 
aircraft are, in general, nonlinear functions of motion variables, their time 
rate of change, and the history of maneuvering. How these unsteady 
aerodynamic forces and moments may be represented becomes uncertain, 
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in particular at high angles of attack. If the response is measured by wind- 
tunnel dynamic testing, questions arise as to how the measured time- 
history data can be analyzed and expressed in a form suitable for flight 
dynamic simulation. For a certain type of nonlinearities produced in a test 
with small- amplitude oscillation, the analysis has been accomplished by 
separating the time-history data into in-phase and out-phase components 
(ref. 14). When large- amplitude forced oscillations are employed in the 
wind-tunnel testing at a large mean angle of attack, the aerodynamic 
phenomena may involve dynamic stall and/or strong vortex flow, with or 
without vortex breakdown. In this case, higher harmonic components in 
the aerodynamic response are expected to exist (ref. 15) and the 
phenomenon of aerodynamic lag may be important. Therefore, a more 
general modeling technique is needed. 

In this research, a numerical method will be developed to analyze the 
nonlinear and time-dependent aerodynamic response to establish the 
generalized indicial function in terms of motion variables and their time 
rates of change. 
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Chapter 2 

Mathematical Development 


2.1 Aerodynamic Modeling 

2.1.1 Historical Review 

For the conventional airplanes, linear relationships are used to 
represent the aerodynamic responses to motion variables (ref. 16). For 
example, the total lift for an airplane in symmetric steady flight can be 
written as 


Cl - C L + C L c^b ( > 

0 ot 

where is the angle of attack measured with respect to wing-body mean 
aerodynamic chord and C Lq is the lift when oc wb =0. C L ^ is the lift-curve 
slope of the whole configuration and a is the angle of attack of the zero lift 
line of the airplane as shown in figure 1. This relationship is good only for 
airplanes flying below stalled a. When airplanes fly at high angles of 
attack, as most modem fighters do, the relationships between the 
aerodynamic responses and motion variables become nonlinear and are 
much more difficult to describe. To deal with the aerodynamic responses 
for helicopters during dynamic stall, Wayne (ref. 17) developed an empirical 
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model to predict the dynamic stall lift. In flight simulation, two common 
ways are used to treat high-angle-of attack aerodynamics. One is using the 
tabulated data (ref. 18) and the other is to use a local linearized model 
which form a piecewise continuous fit of the nonlinear response (ref. 19). 

Based on the functional analysis, Tobak and Schiff (ref. 1) is able to 
develop a fundamental formulation of aerodynamic response for arbitrary 
motion. Their methodology is briefly reviewed in the following. By 
assuming that small step changes of a and qfiAr^ at time x, the incremental 
response aC l at time t is written as 


lim AC L (t,x) 
Act-0 Xa 


= C L [a(0, q«)i t, x] 


lim AC L (t,x) 
A(q(/vJ-0 A(q(/v j 


C L [a(0, q(Oi t, x] 

4 


(2) 


where ^ is a running variable in time over the interval zero to x(see fig. 2), 
f is a reference length and v^ is the free stream velocity. The incremental 
response aC l due to a step change is called the indicial response. If the 
variation of variables a and q over the range 0 to t is replaced by the 
summation of many step changes as shown in figure 3, the summation of 
incremental responses yields an integral form for C L at time t as 
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C L (t) = C L (0) + /;C L [a(0, q(5)i t, t] ^ d t 

•'O • dr 


( 3 ) 


* — /'C L [a(0, q(E); 

v_ J o 4 dx 

To have practical applications, this functional integral form needs to 
be simplified. By assuming that a and q are analytical functions in a 
neighborhood of £=x, variables a and q can be expanded by their Taylor 
series at ^=x. The indicial responses C L , for example can be expressed 
as 


C L.[ a (5)> <1(0; t, x] = C L [t, x; a(x), d(x), , ^ 

q(t), q(x), 1 

By further assuming that only the first two coefficients are needed 
in Taylor expansion of indicial responses, the integral form of eq. (3) 
becomes 


C L (t)=C L (0) + //c l [t,x;a(x),d(x),q(x),q(x)]-^dx 

' dT (5) 

+ — / n ‘ C L [t,x;a(x),d(x),q(x),q(x)]-^dx 
v m Jo q dx 

Eq. (5) is applicable to the study of rapidly varying maneuvers, where 
hysteresis phenomena are known to exist. However, it is difficult to 
implement eq. (5). By introducing the assumption of no hysteresis effect 
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and a slowly varying motion, Tobak and Schiflf neglected the dependence of 
the indicial response on a and q . By further assuming that the indicial 
response is a function of elapsed time t-x instead of t and x separately, a 
much simplified expression of eq. (3) can be written as 


C L (t)=C L (0) + rC L> [t-T;a(T),q(T)]^^d-r 


dx 


-f l C L [t-T;a(T),q(T)]-^^dT 

JO a dt 


(6) 


Although the form of eq. (6) represents a great simplification over that of 
eq. (3), the equation still includes the full linear form as a special case. 

Jenkins (ref. 20) applied a local Taylor expansion to indicial response 
C L and used that Taylor expansion form to fit numerical indicial responses 

a 

calculated from a program called NLWAKE. By substituting C L into eq. 
(6), Jenkins was able to predict the oscillating motion for airfoil at low 
frequencies. 


2.1.2 Current Development 

In the current research, the hysteresis effect is included and the 
assumption of low frequencies will be removed. Therefore, a form between 
eq. (5) and eq. (6) is written as 
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(7) 


C L (t)=C L (0) + / o t C L< [t-T; «(x), a(x), q(x), <*(x)] ^^dx 

+— /‘C, [t-x; a(x), d(x), q(x), 4(x)] ^^dx 

V_ J 0 1 Qx 

In wind-tunnel testing, q is the same as a . Since the method 
developed in this study will be used to analyze wind tunnel data, Ocwill be 
used instead of q in the following investigation and the investigation will 
be focused on lift force. The effect of &(i.e. q ) is included in the virtual 
mass effect (noncirculatory response). As explained in reference 2 for 2-D 
incompressible flow, the noncirculatory response is identical in every 
harmonic motion and therefore, is indepedent of the time history of motion. 
The same concept is adopted in the present study. Then eq. (7) is rewritten 
as 

C L (t)=C L (0)+ noncirculatory response 

*/„' c i, [t - x; “<*>• *< t >] ( 8 ) 

“M. “Ml 

The main objective in the present investigation is to find a suitable form for 
the integrand of eq. (8). Then the time response C L (t) can be calculated 
through the integration by substituting the suitable form of C L and . 
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In the linear theory (refs. 2 and 3) , the aerodynamic response could 
be separated into a product of an amplitude function and a phase function 
in harmonic motion. The amplitude function depends on motion variables 
and their time rate of change. On the other hand, the phase function is a 
function of frequency and accounts for any phase lag between the response 
and the excitation. In a two-dimensional linear theory, the phase function 
is given by Theodorsen’s circulation function (refs. 2 and 3). After response 
has been obtained at different frequencies with the same amplitude in 
harmonic oscillation, the phase function can be determined numerically. 
After use of reciprocal relations (ref. 21), the indicial function can be 
defined by numerical means. This approach has been used for numerical 
determination of indicial lift for plunging airfoil in ref. 5 and for plunging 
wings in ref. 22. 

The method for the linear theory will be generalized as follows. 
Instead of assuming that the aerodynamic response is a product of an 
amplitude function and a phase function as it is in the linear theory, it is 
taken to be a sum of the products of amplitude functions and phase 
functions in harmonic motion; i.e., 

C L =C 0 +E(amplitude function) j * (phase function)^ 
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In the linear theory, j equals 1 in the equation. To determine what 
the forms of the amplitude functions and the phase functions are, the 
aerodynamic response due to harmonic oscillation is assumed to be of the 
form 

C L = F 0 +F(a,dt )a + G(a,d )d (9) 

and it is defined that 

a i = ««. + Oo cos(kt’) 
a = ocq cos(kt’) 

A = (-Ogk) sin(kt’) 

where k is the reduced frequency, t’ is the nondimensionalized time, o^, is 
the mean angle of attack and Oq is the amplitude of angle of attack. To find 
the constant F 0 and functions F and G as functions of a(t) and a (t), a 
functional analysis is needed. However, the following method, "successive 
Fourier analysis," represents a practical way to accomplish the task. The 
first step is to Fourier-analyze the response over one period. For simplicity, 
a Fourier series with three terms will be used to explain the procedure of 
the modeling. Then 
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C L = Afl + A x cos0 + A 2 cos20 + A 3 cos30 
+ B! sin0 + B 2 sin20 + B 3 sin30 


( 10 ) 


The second step is to split the result into the form of eq. (9) with F(a,a ) 
and G(a,a ) being Fourier-analyzed again. The result after "successive 
Fourier Analysis" becomes 

C L = Ao + { CC[0,0] + CC[l,0]a + CC[2,0]a 2 

+ DC[0,l]d + DC[l,l]aa + CC[0,2]d 2 }a 
+ { CS[0,0] + CS[l,0]a + CS[2,0]a 2 

+ DS[0,l]d + DS[l,l]ad + CS[0,2]d 2 }d (11) 


The detailed procedure of "successive Fourier analysis" is shown in 
Appendix 1. By collecting the same order terms of a, a and their products 
together, the result of C L becomes 

C L = Ao + { CC[0,0]oc + CS[0,0]ot } 

+ { CC[l,0]a 2 + DC[0,l]a<A + CS[l,0]ad 

+ DS[0,1] cx 2 } 

+ { CC[2,0]a 3 + DC[l,l]a 2 d 
+ CC[0,2]acx 2 + CS[2,0]a 2 d 

+ DS[l,l]acx 2 + CS[0,2]cx 3 } (12) 
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Equation (12) is a useful form for determining stability derivatives based 
on forced oscillation tests. However, it can be seen that for each different 
frequency k with the same amplitude, there will be different response C L 
and different coefficients CC, DC, CS and DS. To have practical 
applications, a general representation of these coefficients as a function of 
reduced frequency at a constant amplitude is needed. For applications to 
a maneuvering aircraft, the following representation of aerodynamic 
response based on the generalized indicial lift concept is more convenient. 

It is recalled that in the classical airfoil theory the circulatory lift is 
written as the product of Theodorsen’s circulation function (representing 
the phase lag) and the quasi-steady lift (representing the amplitude 
function). In the present nonlinear theory, the same form will also be 
adopted. From the classical potential theory (ref. 23), it has been found 
that Pad6 approximants provide an accurate approximation of the 
theoretical phase function. Therefore, Pade approximants will be used for 
the present model as phase functions to represent coefficients CC, DC, CS 
and DS. Question arises as to what the amplitude function would be. 
Since the quasi-steady lift (lift without lag) represents the amplitude 
function, one way to find the amplitude function is to extrapolate 
coefficients CC, CS, DC and DS from very small k values to k=0. However, 
difficulty may arise in extrapolation as illustrated in the following. 
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From the classical airfoil theory (ref. 2), the lift coefficient for a flat 
plate oscillating with one radian amplitude and with respect to midchord 
in incompressible flow can be written in complex form as 

c g = 2 k [0.5 tit + C(k) (a + 0.5 tic )] = c^ e ,kt 

and 

c*"= 2 k [0.5 ik + C(k) (1 + 0.5 ik)] 

= 2k [0.5 ik + ( F(k) + iG(k) ) (1 + 0.5 ik)] (13) 

where the first term inside the brackets is the virtual mass effect. C(k) in 
the second term is the Theodorsen’s function for the phase lag and the term 
(1 + 0.5 ik) represents the amplitude function or the quasi-steady lift. 
Therefore, the response Cg can be calculated by taking C(k) value from a 
table (ref. 24) and substituting k and C(k) into eq. (13). Suppose the 
responses are known for each k value but not the amplitude function. Then 
the amplitude function needs to be extrapolated from small k values. Eq. 
(13) is rewritten for the unknown amplitude function as 

Cg = 2 k [0.5 ik + C(k) (1 + Hi ik)] (14) 


14 



Separating eq. (14) into the real and imaginary parts, it is obtained as 


= 2tc (F(k) - H x G(k) k) 


(15a) 


and 

= 2k (0.5 k + G(k) + H t F(k) k) (15b) 

For k=0.01, the exact values of the phase function (ref. 24) are 
F(k)=0. 9824215 and G(k)=-0. 0456521. The response Cg is calculated from 
a very accurate 2-D UQVLM program and <^=-0.2245678. Then, is 
obtained by dividing the small difference 

^ 2 k - (0.5k + G) = 0.0049110 

by a small value of k and 1^=0.49989. In application to a general problem, 
G, F and Hj are all unknowns. Therefore, numerical errors will be 
exaggerated in predicting H : and the result is unpredictable. For nonlinear 
cases, the phase functions are also unknowns and need to be solved. 
Therefore, a small error in assuming the phase function will make the 
prediction of quasi-steady response inaccurate. The numerical solution in 
a direct optimization method would be difficult to converge as has been 
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experienced. To avoid this kind of difficulty, a form for the amplitude 
function is modeled in the present research and to fit the aerodynamic 
responses at given k values by a one dimensional gradient method. For 
this purpose, eq. (9) (or the experimental oscillatory results) is rewritten in 
a complex form, as follows: 


C L = Ao + (A, - iBj) e ikt ’ + (A* - iB 2 ) e i2kt ' 

+ (A, - iB 3 ) e ,3kt ’ ( 16 ) 

It should be kept in mind that only the real part of the response has 
a physical meaning. The reason to put in the complex form is to benefit 
from the great mathematical convenience of the e ikt ' notation. If a is 
rewritten as 


a = clq e ikt ‘ 


and 


a = (icxok) e ikt ' 
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then eqs. (12) and (16) and the classical airfoil theory suggest that the 
response could be put in the following form involving the products of 
amplitude functions and phase functions as 

C L - C 0 (k) 

+ E u d + Ejjfi + C, * (H n o + H^d ) * (1 - PD,) 

+ E 12 d 2 + E 22 & 2 + q * (H, 2 a 2 + Eqad + H^d 2 ) 

* (1 - pd 2 ) 

+ E 13 d 3 + E 23 a 3 + C 3 * (H 13 a 3 + H 23 a 2 a + FL^ad 2 + H 43 a 3 ) 

* (1 - PDj) 

where PD is a Pad6 approximant with order 2 and is defined as 


p D p n ^) 2 + p 2i Cik) 

J P 3J (ik) 2 H- (ik) + P 4j 

E n a, + E 21 G.J etc. are the virtual-mass effect and account for the 
noncirculatory lift (ref. 2). The variables a i and a , are defined as 


a j = ik aj e 1 


ijkf 


and 
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tt j = -k 2 Oq e ijkt ’ 


to be consistent with higher order terms. When j=l in the above equations, 
a j = a and & j = ft . In addition, H 21 , H 22 , H^, etc., are related to the pitch- 
rate effect. It should be noted that those terms inside the parentheses 
following C lf C 2 , C 3 , such as (H n a + H 21 a ), represent the quasi-steady 
response and (1 - PDp represents the unsteady aerodynamic lag in 
response. Therefore, the present assumed form for aerodynamic modeling 
encompasses the classical linear theory. 

Cj are the reference values used to normalize the lift given by - i 
Bj in the least squared-error method, j is the index consistent with the 
exponent of the exponential term in eq. (16). For example if the j’s term in 
eq. (17) represents the coefficient of e lkt , then j is 1. If the j’s term in eq. 
(17) represents the coefficient of e l2kt then j is 2, etc. The first term, C 0 (k), 
in eq. (17) is a constant term, supposedly a function of frequency. From 
available experimental data (ref. 25), it is found that an averaged constant 
can be used to represent C 0 (k) term as shown in Figure 4 for a delta wing. 
The unknown coefficients Py, P^, P 3j and P^ are calculated from the least 
squared-error method. E n , E 21 , H u , H 12 , etc., will be obtained separately 
by minimizing the sum of squares of errors. This is equivalent to a two- 
level optimization method to determine the unknowns in eq. (17). That is, 
E, H, etc., are assumed first. Then P y , etc., are determined by minimizing 
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the sum of squared errors. The values of E llf etc., are varied next so 
that the sum of squared errors is minimized. It was found that this 
approach is more effective in determining a global minimum solution for 
the unknowns than a straightforward optimization (one level) method 
because of nonlinearity in the unknowns in the optimization problem. It 
should be noted that in the literature the phase function has been typically 
determined by the response to plunging motions. Therefore, those terms 
associated with a in eq. (17) do not appear. This would very much simplify 
the mathematics of determining the Pade approximants. The details of the 
present method are discussed in the following. 

2.2 Least-Sauare Method 

By choosing proper values of E n , H u , H 12 , etc., in eq. (17), the 
corresponding - i Bj term in eq. (16) is then divided by the amplitude 
function. The result will appear as 


i _ Aj - iBj - E tj ik - Eg (-k 2 ) 

(amplitude function), 

J (18) 

P lf (ik) 2 + P 2j (ik) 

P 3j (ik) 2 + (ik) + P 4j 


If both sides of eq. (18) are multiplied by the denominator of the Pad6 
approximant and separated into real and imaginary parts, then 
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Re = P y k 2 - P^k 2 + Py Vj - Wj k = 0 


(19a) 


and 


Im S P^k + P 3j Wjk 2 - PyWj - Vjk = 0 


(19b) 


The sum of squared errors is defined as 

Err ■ E Re(kj) 2 + E Im(kj) 2 (20) 

By equating the first derivatives of squared errors (eq. 20) with respect to 
variables P y> P^, P 3j and P 4j to zero, the unknown coefficients P, , P, P,. 
and P 4j can be determined by 

Eq 0 -EVjkf EVjlq 

0 Elq EWjkf -EWjkj 

"EVjkj 4 EWjkf ECV^.wjV) -WK+WX) 

Ev.k 2 -EWjkj E(vf kf+wft 2 ) Ecvf+Wj 2 ) 

where i varies over the range of input frequencies, and the mode subscript 
j on V and W has been omitted. 


2j 

P 3j 


EWjkf 


0 

0 


( 21 ) 
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2.3 Gradient Method 


After the unknown coefficients Py, P^, P 3j and P 4j have been found, 
a one-dimensional gradient method is used to find E and H values which 
will make the sum of the squared errors minimum The E or H value is 
perturbed first by a small amount AE or AH to find the gradient of the s um 
of squared errors. If the gradient tends to reduce the error, then the E or 
H value is perturbed further until several iterations has been reached (it 
is set to be 5 iterations in the current program). After that, the same 
procedure is applied to other E or H. Then the whole procedure is repeated 
again until several iterations has been reached (it is set to be 10 in the 
current program). 

Finally, the response of C L is written as 
C. = C 

L *ve 

+ E n d + E^a + C, * ( H n a + H^d) * (1 - PD,) 

+ E 12 d 2 + + C 2 * ( H, 2 a 2 + H^ad + H^d 2 ) * (1 - PD 2 ) (22) 

+ E 13«3 + E 23®3 

+ C 3 * ( H 13« 3 + + H3 3 od 2 + H 43 d 3 ) * (1 - PDj) 

It is easily seen that each term in the above equation is a product of an 
amplitude function and a phase function. The procedure to put oscillating 
response data into the form of eq. (22) is summarized in the next section. 
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2.4 Summary of Numerical Procedure 


Step 1. Unsteady-response analysis: 

Use Fourier analysis to analyze the harmonic motion responses 
for different frequencies over one period. For each frequency, the responses 
should be in the same form as in eq. (16). 

Step 2. Constant-term analysis: 

From step 1, if constant terms appear in the Fourier analysis 
then C ave is calculated from each constant term due to different frequencies 
as 


N 

E Ao^n) (24) 

n* 1 


where N is the number of frequencies used in step 2 and are the 

constant terms due to different frequencies k n in Fourier analysis. 

Step 3. Amplitude-phase identification: 

In this step, the coefficients Aj and Bj calculated from step 1 are put 
into the form of a product of amplitude functions and phase functions as in 
eq. (22). The procedure is as follows: 

3-a Set the reference values Cj 
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3-b Set the initial guess for E or H values. 

3-c Use the least-square method to find unknown coefficients 
^ij» f*3j an d P 4 j. 

3-d Use the gradient method to find better E or H values. 

Repeat steps 3-b to 3-d until the sum of square errors reached 
minimum or the iteration limit is reached. 

Although three-term Fourier series are used in the above, the 
procedure is applicable to any number of Fourier terms. 

2.5 Indicial Formulation 

In linear theory, the reciprocal relations (or Fourier summation) (ref. 
21) has been used to calculate the indicial response. However, in nonlinear 
theory those relations can not be applied. As Tobak (ref. 7) mentioned in 
his paper, the aerodynamic response due to a step change should reach 
steady-state value asymptotically at subsonic speeds. In linear theory, 
these asymptotic relations are represented by exponential functions (refs. 
2 and 5), and these exponential functions are calculated through the phase 
function. Therefore, the phase function in the current nonlinear modeling 
will be converted into exponential function in time domain but the 
amplitude function is kept unchanged. If eq. (21) is rewritten for m-terms 
Fourier series as 
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+ E u d + Cj * (H n a + H^d) * (1 - FD^ 

+ E 12 d 2 + E 2Z a 2 * C 2 * (H 12 a 2 + H^ad + H3 2 d 2 ) * (1-PD 2 ) 
+ E l3 d 3 + Ejjdj + C 3 * (H 13 a 3 + H^d + H»«d 2 + H 43 d 3 ) 
* (l-PDj) 

+ 

c + 

w «ve 

m 

+ E Eyffj + ^2j®j 
J-i 

+ (amplitude function) j * (phase function ) } 


( 25 ) 


where the phase functions represented by the Pade approximants are 
defined as 


j _ P„ (ik ) 2 + P 2j (ik) = i _ ik aij + ika^, 

P 3J (ik ) 2 + ik + P 4j ik + a 3j ik + a 4j (26 ) 

_ i(jk) a }j i(jk) a 2j 

i(jk) + ja^ + i(jk) + ja 4J 

By applying the Fourier integral to the phase function, the corresponding 
expression in time domain is calculated as 
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(27) 


1 

2* i(jk) 



p ij(ft) 2 + PaOk) 

P3j(ik) 2 + ik + P 4J 


] e^'d (jk) 




- *2i e 


V 


the detailed derivation is in appendix B. Then, the indicial response C L for 
the circulatory lift is written as 


C l_. " c i *(H„a + H^a) * (1 - a n e'* 31 ' - 
+ c 2 * (Hu® 2 + H^ad + Hjjd 2 ) 

* 0 - *k e”” 21 ' - *1, e ^ 1 ') 

* (28) 

m 

- £ (amplitude function), 
j-i J 

* (1 - a y e"*^ - e"^) 

In the current research, five terms Fourier series is used to analyze 
the experimental data. Therefore, the response can be written as 


- Ci * (H n a + H 21 d) * (1 - exp T ) 

+ C 2 * (H 12 a 2 + H 22 ad + H 32 d 2 ) * (1 - exp 2 ) 

+ C 3 * (H 13 a + H 23 a 2 a + H 33 ad 2 +• H^d 3 ) * (l-exp 3 ) 

+ C 4 * (H 14 a 4 + H 24 a 3 d + H 34 a 2 d 2 + H^ad 3 + H,,d 4 ) (29) 

* (1 - exp 4 ) 

+ C 5 * (H 16 a 5 + H^o^d + H 36 a 3 d 2 + H 46 a 2 d 3 + H 56 ad 4 
+ H 65 a 5 ) * (1 - exp 5 ) 
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By differentiating with a and a , are obtained 

as follows: 


- c, • H n . (1 - op,) 

+ C 2 * (2H l2 a + H^d) * (1 - exp^ 

+ C 3 * (3H 13 a 2 + 2H2 3 ad + H^d 2 ) 

* (1 - expj) 

+ C 4 * (4H 14 a 3 + 3H 24 a 2 d + 2H 34 ad 2 + H^d 3 ) 

* (1 - exp^ 

+ C 5 * (5H 15 a 4 + 4H2 5 a 3 d + 3H 35 a 2 d 2 + 2H 45 ad 3 

* (1 - expj) 


(30) 


H55* 4 ) 


(C Lw J* - Cj • H„ * (1 - exp,) 

+ C 2 * (H^a + 2H 32 d) * (1 - exp 2 ) 

+ C 3 * (H^a 2 + 2H3 3 ad + 3H 43 d 2 ) * (1 - expj) 

♦ C 4 . (H^a 3 * 2H 34 a 2 a - 3H„ad 3 * U^d 3 ) (31) 

* (1 - expj 

+ C 5 * (H^a 4 + 2H3 5 a 3 d + 3H 45 a 2 a 2 + 4H 55 ad 3 + SH^d 5 ) 

* (1 - exp^ 

By substituding these derivatives into equation (8), the integrand can 
be determined. The initial condition C L (0) is calculated by setting the 
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variable x to zero in the indicial response C L . Since the noncirculatory 

indidal J 

lift is identical in every harmonic motion, it must be excluded from the 
integral. Then the generalized response for arbitrary motion is obtained by 
substituted the noncirculatory lift and the derivatives of the indicial 
response to eq. (8). Therefore, the final form of time integration for 
arbitrary motion is written as 


C L (t/) = “WUo + C .vc + £ (Eijttj + 

£ /V d (amplitude function). 

to MO d . 

♦ -Lzf d (amplitM<ie fupcti ° n >' ,n „ e -V-' -> . e-V'-'X d*(T) 
V^-1 J 0 d a U *j e V } ~ZT 


dx 


( 32 ) 


2.6 Arbitrary motions 

In the current research, the indicial responses are derived from data 
of harmonic motions. Therefore, to have a correct representation, an 
arbitrary motion must be represented locally by a cosine function. Since 
the values of a x and a are usually given at a certain time, they can be 
described by the cosine and sine function as 
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= «m + Oo cos(<(> + kt’) = a m + a 


( 33 ) 

a = -a,, k sin(<|) + kt’) 

By defining a mean angle-of-attackCo^) and amplitude^), the above 
equations can be rewritten as 


F : = a - ocq cos(<|) + kt’) = 0 


F 2 = a + oto k sin(<]) + kt’) = 0 


(34) 


the unknowns, the reduced frequency k and phase function angle <() at a 
given time t’, can be solved by Newton’s method(see Appendix C for a 
detailed derivation) 
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Chapter 3 

Results and Discussion 


3.1 Linear Results 

Because appropriate high-alpha experimental data to apply the 
present modeling method are limited, the present method will first be 
tested with linear theoretical results. Several cases in the two-dimensional 
and three-dimensional linear flow have been studied to verify the proposed 
method of aerodynamic modeling. 

3.1.1 Two-Dimensional Flow 

The first case studied is a 2-D flat plate oscillating in the 
incompressible flow. The amplitude is 57.3 degree (one radian) in angle of 
attack for the airfoil. Therefore it oscillates from 57.3 degree of angle of 
attack to -57.3 degree of angle of attack then back to 57.3 degree of angle 
of attack for one cycle with respect to midchord, i.e. 

a! = 0.0 + 1.0 cos(kt’) (in radian) 

The steady lift is already known from the linear theory (ref. 2) as 
27tcc, and the oscillating complex lift is taken from a 2-D unsteady QVLM 
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program (ref. 26) as input data for the current model and for comparison. 
Through numerical experimentation, it is found that only six frequencies 
are needed to have accurate results. Through the modeling procedure as 
summarized in section 2.4, the lift can be written as 

c l = ® + E 21 ix + 2 k ( H n a + H 21 a ) (1.0 - PDj) 

and the values of coefficients E n , E 21 , H n , H 21 and P n (for pade approximant 
PD X ) are listed in Table 1. The results for the lift coefficients are plotted 
in Figure 5 for different numbers of frequencies used. Compared with the 
aerodynamic responses by the 2-D UQVLM program, the numerical results 
from modeling show excellent agreement. 

Two Mach numbers, 0.2 and 0.4, are chosen in the 2-D compressible 
flow to verify the current model. A two-dimensional unsteady QVLM (ref. 
26) program is again used to calculate the complex lift as input data for the 
current model and for comparison. The same frequencies as in the 
incompressible flow are used as input also. The results for coefficients E u , 
^ 2 i> H u , H 21 and P u are listed in Table 1. The aerodynamic responses c^ 
calculated by the model are plotted in Figures 6 and 7 to compare the 
results from 2-D unsteady QVLM. These figures show that the numerical 
results by modeling are very accurate. 
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3.1.2 Three-Dimensional Flow 

The same Mach numbers (includes 0) in the 2-D flow are used in 3-D 
attached flow to verify the current model. The geometry is a 70 degree 
delta wing which oscillate from zero degree angle of attack to twenty degree 
angle of attack with respect to mid root chord, i.e. 

aj = 0.1745329 + 0.1745329 cos kt’ (in radian) 

This means that the mean angle of attack is ten degree (0.1745329 
radian) and the amplitude of the oscillation is ten degree (0.1745329 
radian). The input aerodynamic responses are calculated from a 3-D 
unsteady QVLM program (ref. 27). In the program, the total lift is the 
sum of steady lift at the mean angle-of-attack plus unsteady lift. Since the 
steady lift is the same for every term, only the unsteady lift is used in the 
modeling and for comparison. Through numerical experimentation, it is 
found that the responses at low frequencies do not change significantly, 
which results in inaccurate modeling. To have accurate approximation, 
high frequency responses are needed. Seven reduced frequencies (k = 0.01, 
0.1, 0.2, 0.6, 1.0, 2.0, 2.5) are used as input data in the 3-D attached flow 
cases. The results for the coefficients of the modeling are listed in Table 1. 
The responses C L from modeling are plotted in Figures 8 to 10 to compare 


31 




with the results from 3-D unsteady QVLM program. All of these figures 
show very good agreement. 

3.2 Nonlinear Results 

The experimental data (ref. 25) for a 70-degree delta wing in pitching 
oscillation is used to validate the current aerodynamic model. The angle 
of attack which describes the pitching motion with respect to the 57% root 
chord is given as 

cti = 27.5 + 27.5 cos(180 + kt’) (in degree) 

which means the delta wing oscillates from zero degree angle of attack to 
55-degree angle of attack then back to zero degree angle of attack for one 
cycle. The reduced frequency k is nondimensionalized based on wing’s root 
chord and the pitching moment is measured with respect to quarter root 
chord. Five sets of data corresponding to five different frequencies are 
available and they will be used as the input data to calculate the 
coefficients for the current aerodynamic model. Five terms in the Fourier 
series are employed for the calculation. The calculated coefficients for the 
current aerodynamic model are listed in Table 2 to Table 5 for C L , C D , C m 
and C^. 
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The lift coefficients obtained from the aerodynamic modelCeq. 29) are 
compared with the original test data in Fig. 11 with good agreement. 
Experessions for C D , C m and C N are obtained with the same procedures as 
those used for C L . The modeled harmonic results are compared with data 
in Figs. 12 to 14 respectively. Again, the good agreement indicates that the 
present aerodynamic model is accurate in representing the experimental 
harmonic data. 

3.2.1 Indicial Formulation 

Note that eq. (8) is valid for arbitrary motion. To check its range of 
validity in the nonlinear theory, extensive study has been conducted. 
Because the indicial responses consist of oscillatory responses, two 
oscillatory lift cases in the last section are used to verify the indicial 
integration first. That is, by assuming oscillatory motion in eq. (8), the 
time-integrated lift response should agree with the forced-oscillation data. 
Since the angle of attack is set to be a complex number (cos(kt’)+i sin(kt’)) 
in oscillating cases, only the real part of the integrated lift is taken. The 
lift by integrating eq. (8) for a 70-deg. oscillating delta wing with 
frequencies k=0.098 and k=0.165 are plotted in Figure 15. Compared with 
the lift from aerodynamic modeling(eq. 25), the integrated lift shows good 
agreement. Which confirms that at least the current methodology of 
indicial integration works for harmonic oscillations. 
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The experimental data(ref. 26) available for more general motion is 
for harmonic ramp motion. That is, the 70-deg. delta wing model undergoes 
a harmonic pitching motion(with respect to the 57% root chord) to a certain 
angle-of-attack and then stops there. To verify the aerodynamic models 
further, these harmonic ramp responses are applied. Note that the 
experimental harmonic ramp motion is described by the following pitching 
motion until reaching a certain angle-of-attack. i.e. 

= 27.5 + 27.5 cos(180 + kt’) (in degree) 

The results calculated from indicial integration for C L in ramp motion 
up to 55 deg. angle-of-attack at two different frequencies are plotted in 
figure 16. Compared with experimental data for the 70-deg. delta wing(ref. 
28), the C L responses show good agreement for frequencies 0.0926 and 
0.0714. The results for C L in ramp motion up to a=35,45 and 55 deg at a 
frequency equal to 0.0714 are ploted in figure 17. It is seen that the 
present aerodynamic model is fairly accurate if the harmonic ramp motion 
is from ct = 0 to 55 deg. However, the final C L is overpredicted if the ramp 
motion stops at an a less than 55 deg, even though the peak C L is still well 
predicted. A possible reason for this is that the harmonic data based on a 
= 0-55 deg. contain dynamic effect on vortex-breakdown characteristics at 
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a < 55 deg. The results for C L at a steady a = 35 or 45 deg. should be 
higher than the static values. 

The corresponding drag and pitching moment coefficients at reduced 
frequencies 0.0926 and 0.0714 are presented in figure 18 and figure 19. For 
the C m responses, the result is well predicted except at small time. 
However, the drag coefficient is not as well predicted in ramp motions as 
in harmonic motions(see Fig. 12). The same descrepancy in predicted drag 
coefficients is also found in figure 20 for a harmonic ramp motion with a = 
0 to 35 and 45 deg at a frequency eaual to 0.0714. It is not known whether 
this is caused by differences in the test models and test Reynolds numbers. 
The test model for the harmonic motions(ref. 25) has two-sided chamfered 
leading edges with a thickness of 0.5 inch at a Reynolds number of 1.64x10® 
based on the root chord. The model for the ramp motions(ref. 28) is 
chamfered only on the lower surface of the leading edge and has a thickness 
of 0.25 inch, and tested at a Reynolds number of 1.54x10®. The effect of 
Reynolds number on the static lift coefficient for two different delta wing 
models is taken from ref. 29 and reploted in figure 21. As shown in the 
figure, the effect of Reynolds number for the one-side chamfered delta wing 
appears to be confined to angle of attack above 28 deg and much less than 
for the two-side chamfered delta wing. The effect of Reynolds number on 
dynamic lift is reploted as figure 22. As the Reynolds number is increased, 
the maximun lift tends to decrease. The similar effect of Reynolds n um ber 
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to other dynamic responses of the two-sided chamfered delta wing can be 
found in ref. 25. Therefore, using the harmonic motion data of one model 
to predict the ramp motion of another model may cause discrepancy in drag 
coefficients. 

To illustrate the application of the present aerodynamic modeKeq. 30) 
for arbitrary motions, a linear ramp motion is assumed in the integration. 
As described in section 2.6, the linear ramp motion is represented by an 
equivant local harmonic motion using a cosine function. However, from 
numerical experimentation it is found out that due to the discontinuity in 
slope(a ), if the mean angle and amplitude are set to have the same value, 
then the calculated responses will have discontinuity there. To avoid the 
discontinuity of the responses, the amplitude must be set larger than the 
mean angle of attack. Three values of the amplitude have been tested for 
the linear ramp motion from a = 0 to 55 deg. and the calculated lift 
responses are plotted in figure 23. It can be seen from the figure that a 2.5 
deg increase in amplitude will smooth out the response. Therefore, in 
calculating the linear ramp response the amplitude will be set 2.5 deg. 
larger than the mean angle-of-attack. The lift coefficient for the linear 
ramp motion for cx = 0 to 55 deg is compared with that in a harmonic ramp 
motion in figure 24. It is seen that the linear ramp motion tends to 
produce higher C L beyond the peak value because it has a higher value in 
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& . Similar results can be found in figure 25 for lift coefficients 
corresponding to the linear ramps for a = 0 to 35 and 45 deg. 

The indicial lift responses consist of harmonic data which are based 
on a = 0 to 55 deg are now used to calculate the lift responseto a linear 
ramp motion up to a = 76 deg. As shown in figure 26, the trend is 
predicted but the magnitude is over-predicted. Therefore, it may be 
concluded that the harmonic data based on ct = 0 to 55 deg. may not be 
used to extrapolate the data to higher angles of attack. 

The last case is to test the mean angle-of-attack effect. That is, using 
the harmonic data(ref. 25) at = 27.5 deg., harmonic responses under 

different mean angles-of-attack are calculated. From numerical 
experimentation, it is found that to have correct trend in the prediction, if 
the mean angle of attack is greater than 27.5 deg, it is set to 27.5 deg. The 
results from indicial integration are plotted in figure 27. Compared with 
experimental data taken in a different tunnel with different 
instrumentation(ref. 31), these numerical results show correct trend. Note 
that the experimental data are taken under Re = 1.64xl0 6 based on the root 
chord(4.166 ft) for a 70-deg delta wing. The Reynolds number is the same 
as that of the harmonic data with 27.5 deg. mean angle-of-attack for the 
aerodynamic model but based on a larger size model and smaller free 
stream velocity. In addition, the larger 70-deg model is oscillating with 
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respect to the 40% mean aerodynamic chord(60% root chord) instead of 57% 
root chord. Therefore, it is difficult to assess the accuracy of the present 
aerodynamic model. 
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Chapter 4 

Conclusions and Recommendations 


4.1 Conclusions 

In the current study, a method based on Fourier analysis is 
developed to analyze the force and moment data obtained in large- 
amplitude forced oscillation tests at high angle-of-attack. Results from 
linear two- and three- dimensional unsteady aerodynamic theories as well 
as test data for a 70-deg delta wing are used to verify the models. It is 
shown that the present modeling method is accurate in producing the 
aerodynamic responses in C L , C D , C m and C N to harmonic motions and the 
ramp- type motions. The model also produces a correct trend for a 70-deg 
delta wing in harmonic motion with different mean angles-of-attack. To 
deal with the linear ramp type motions, a local harmonic motion represen- 
tation is used in the current model. However, the current model failed to 
predict the response to a ramp motion which has larger angles-of-attack 
than that of the harmonic motions which form the aerodynamic model. 
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4.2 Recommendations 


In the current model, the effects of Reynolds number and Mach 
number are not included in the modeling. The main reason is that not 
enough experimental data or theoretical results are available to support the 
study. Therefore, the recommendations for the future work are 

(i) A consistent experimental work or numerical calculation for unsteady 
longitudinal aerodynamics should include Reynolds number and Mach 
number effects for each reduced frequencies in harmonic motions. 

(ii) Modify the current mathematical model to include the Reynolds number 
and Mach number effects based on the results from (i). 
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Figure 4 Comparison of Numerically Approximated Constant 
Terms with Results from Fourier Analysis for 
The lift of an Oscillating 70— Deg. Delta Wing 



Fig. 3 Summation of incremental response 
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Lift Coefficient C, 



Reduced Frequency k 

Figure 5 Unsteady Lift Coefficient For a 2— D Flat Plat 
Pitching about Midchord at M = 0.0 
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Lift Coefficient C, 



Figure 6 Unsteady Lift Coefficient For a 2— D Flat Plat 
Pitching about Midchord at M = 0.2 
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Lift Coefficient C, 



Figure 7 Unsteady Lift Coefficient For a 2-D Flat Plat 
Pitching about Midchord at M = 0.4 
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Lift Coefficient C 



Figure 8 Unsteady Lift Coefficient For a 70— deg Delta Wing 
Pitching about Mid— Root chord at M = 0.0 
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Lift Coefficient C 



Figure 9 Unsteady Lift Coefficient For a 70— deg Delta Wing 
Pitching about Mid-Root chord at M = 0.2 
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Lift Coefficient C 



Figure 10 Unsteady Lift Coefficient For a 70— deg Delta Wing 
Pitching about Mid-Root chord at M = 0.4 
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Lift Coefficient C. Lift Coefficient C 




Angle of Attack 

Figure 1 1 Unsteady Lift Coefficient for a 70— deg Delta Wing 
Pitching about 57% of Root Chord at Low Speed 
and Various Frequencies 
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Figure 11 Continued. 
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Lift Coefficient C. Lift Coefficient C 



Angle of Attack 
Figure 1 1 Concluded. 


54 



Drag Coefficient C Q Drag Coefficient C 



Figure 12 Unsteady Drag Coefficient for a 70-deg Delta Wing 
Pitching about 57% of Root Chord at Low Speed 
and Various Frequencies 
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Figure 12 Continued. 
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Drag Coefficient C Q Drag Coefficient C 



Angle of Attack 
Figure 12 Concluded. 
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Pitching Moment Coefficient C Pitching Moment Coefficient C 




Angle of Attack 

Figure 13 Unsteady Pitching Moment Coefficient for a 70— deg 
Delta Wing Pitching about 57% of Root Chord at 
Low Speed and Various Frequencies. Moment 
Center at 25% of Root Chord 
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Figure 13 Continued. 
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Pitching Moment Coefficient C Pitching Moment Coefficient C 




Angle of Attack 
Figure 13 Concluded. 
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Normal Force Coefficient C N Normal Force Coefficient C 




Figure 


14 


Angle of Attack 

Unsteady Normal Force Coefficient for a 70-deg Delta 
Wing Pitching about 57% of Root Chord at Low Speed 
and Various Frequencies 
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Figure 14 Continued. 
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Figure 14 Concluded. 
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Nondimensional Time t* 

Figure 15 Unsteady Lift Coefficients from Numerical Modeling 
and Indicial Time Integration for a 70— deg Delta 
Wing in Harmonic Pitching Oscillation about 
57% of Root Chord at Low Speed 
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Nondimensional Time t' 

Figure 16 Unsteady Lift Coefficients from Indicial Lift Model 
and Experiment for a 70— deg Delta Wing in 
Harmonic Ramp Motion at Low Speed 
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Figure 17 Unsteady Lift Coefficients from Indicial Lift Model 
and Experiment for a 70-deg Delta Wing in 
Harmonic Ramp Motion at Low Speed 
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Figure 18 Unsteady Drag Coefficients from Indicial Drag Model 
and Experiment for a 70-deg Delta Wing in 
Harmonic Ramp Motion at Low Speed 
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Figure 18 Concluded. 
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Figure 19 Unsteady Pitching Moment Coefficients from Indicial 
Pitching Model and Experiment for a 70-deg Delta 
Wing in Harmonic Ramp Motion at Low Speed. 
Moment Center at 25% of Root Chord 
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Figure 19 Concluded. 
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Figure 20 Unsteady Drag Coefficients from Indicial Drag Model 
and Experiment for a 70-deg Delta Wing in 
Harmonic Ramp Motion at Low Speed 
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Figure 22 Effect of Renolds Number on Dynamic Lift Coefficient 



Angle of Attack 

Figure 21 The Effect of Reynolds Number on The Static Lift 
Coefficient for Two Different Delta Wing Models 
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Figure 23 Unsteady Lift Coefficients from Indicial Lift Model 
for a 70— deg Delta Wing in Linear Ramp Motion 
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Figure 24 Unsteady Lift Coefficients from Indicial Lift Model 
for a 70— deg Delta Wing in Both Linear and 
Harmonic Ramp Motion 
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Harmonic Ramp Motion at Low Speed 


77 


Angle of Attack (Deg.) 


Normal Force Coefficient C 



78 



Normal Force Coefficient C N Normal Force Coefficient C 




Angle of Attack 
Figure 27 Concluded. 
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Table 1: Results for Linear Flow 
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Table 2: Coefficients for C L in Harmonic Motion 
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Table 3: Coefficients for C D in Harmonic Motion 
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Table 4: Coefficients for C m in Harmonic Motion 
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Appendix A 

Successive Fourier Analysis 


The first step of "successive Fourier analysis" is to Fourier-analyze 
the response over one period. For simplicity, a Fourier series with three 
terms will be used to explain the procedure of the modeling. Then 

= Ao + Aj cosQ + A 2 cos29 + A 3 cos30 

+ B 1 sin0 + B 2 sin20 + B 3 sin30 (Al) 


where 


c ^ d9 

A n = 1 P 2 * C L cos(n0) d0 

K 

C L sin(n0) d0 
n=l,2,3 and 0=kt' 


(A2) 


Once the coefficients of Ao, A„ and B n have been found, the next step 


is to split the coefficients into two groups by using the following formulas, 


A.1 



cos n0 = C(n,0) cos n 0 - C(n,2) cos n - 2 0 sin 2 © 

+ C(n,4) cos n ' 4 © sin 4 0 + 

(A3) 

sin n0 = C(n,l) cos 0 ' 1 © sin0 - C(n,3) cos“- 3 0 sin 3 0 
- C(n,5) cos 0 ' 8 © sin 5 0 + 


where 

C(n,m) = 1 and n!=l*2*3*4******n 

(n-m)! m! 

Therefore, the response of C L becomes 

= Aq + Aj cos0 + A 2 [cos 2 0 - sin 2 0] 

+ A 3 [cos 3 0 - 3cos0 sin 2 0] 

+ sin0 + B 2 [2cos0 sin0] 

+ B 3 [3cos 2 0 sin© - sin 3 ©] 

= Aq + [A t + A 2 cos© + A^cos 2 © - 3sin 2 ©)] cos© 

+ [Bl + 2B 2 cos© + B 3 [3cos 2 © - sin 2 ©] sin© 

= Ao + F(cos© , sin©) cos© + G(cos© , sin©) sin© 

Perform the Fourier analysis again for functions F(cos© , sin©) and 
G(cos© , sin©) by using Fourier series with the same terms as in the first 


A. 2 



step. Then 


F(cos9,sin0) = F 0 + FA, cosG + FAj cos20 + FA 3 cos30 

+ FB, sin0 + FB 2 sin20 + FB 3 sin30 
G(cos9,sin0) = GO + GA, cos0 + GA;, cos20 + GA 3 cos30 

+ GB, sin0 + GB 2 sin20 + GB 3 sin30 

where 


II 

o 

fc 

1 
2 K 

jT Fde 

II 

o 

o 

1 

~2k 

f" g d 6 

FA n = 

h 

fFde 

GA n = 


^ G d 8 

0 

FB n = 

tJ 

P 2 ” F d 0 
0 

GB n = 

U 

G d 0 

0 


n = 1,2,3 


Using eq. (A3) again, then 

F(cos0,sin0) = F 0 + FA, cos0 + FA* [cos 2 0 - sin 2 0] 

+ FAa [cos 3 0 - 3 cos0 sin 2 0] 

+ FB, sin0 + FB a [2cos9 sin0] 
+ FB 3 [3cos 2 0 sin0 - sin 3 0] 

G(cos0,sin0) = G 0 + GA, cos0 + GA,, [cos 2 0 - sin 2 0] 
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+ GA 3 [cos 3 0 - 3 cos0 sin 2 0] 

+ GBj sin0 + GB 2 [2 cos 0 sinO] 

+ GBg [3cos 2 0 sin0 - sin 3 0J 

Therefore 

Cl = Aq + {F 0 + FA! cos0 + FAg [cos 2 0 - sin 2 0] 

+ FA 3 [cos 3 0 - 3 cos0 sin 2 0] 

+ FB, sin0 + FB 2 [2 cos 0 sinO] 

+ FBg [3cos 2 0 sinO - sin 3 0]} cos0 
+ {G 0 + GAj cos0 + GA 2 [cos 2 0 - sin 2 0J 
+ GAg [cos 3 0 - 3 cos0 sin 2 0] 

+ GB! sinO + GB 2 [2 cos 0 sin0] 

+ GBg [3cos 2 0 sin0 - sin 3 0]} sinO 

All the terms associated with cos n 0 on the right hand side of the 
above equation are divided by (a^) D and the terms associated with sin n 0 are 
divided by (-k<Xo) n . After rearrangement, the response of C L becomes 

C L = Ao + {CC[0,0] + CC[1,0] a + CC[2,0] a 2 + CC[3,0]a 3 
+ DC[0,1] cx + DC[1,1] act + DC[2,1] a 2 d 
+ CC[0,2] <x 2 + CC[1,2] adx 2 + DC[0,3] <x 3 } a 
+ {CS[0,0] + CS[1,0] a + CS[2,0] a 2 + CS[3,0]a 3 
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+ DS[0,1] a 


+ DS[1,1] act + DS[2,1] a 2 fc 


+ CS[0,2] Cc 2 + CS[1,2] act 2 + DS[0,3] dt 3 } a (A4) 


where 


CC[n,m] 

CS[n,m] 


FA. 


[a 0 n (-ka 0 m )] 
[<*»" (-ka”)]' 


DC[n,m] 

DS[n,m] 


FB. 


[a" (-ka 0 m )j 


GB„ 


[a 0 n (-ka 0 m )] 


and the coefficients CC[n,m], DC[n,m], CS[n,m] and DS[n,m] are zeros for 
n+m i 3. Comparing with eq. (Al), it is obtained that 


F 0 - Aq 

F(a, Ct) = CC[0,0] + CC[1,0] a + CC[2,0] a 2 

+ DC[0,1] a + DC[1,1] aa + CC[0,2] Ct 2 
G(a, Ct) = CS[0,0] + CS[1,0] a + CS[2,0] a 2 

+ DS[0,1] cx + DS[1,1] acx + CS[0,2] cx 2 
Finally, collecting the same order terms together, then 

Gl = Ao 

+{ CC[0,0]a + CS[0,0]& } 
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+{ CC[l,0]a 2 + DC[0,l]aa 


+CS[l,0]otd: + DS[0,1] <jfc 2 } 
+{ CC[2,0]a 3 +DC[l,l]a 2 c* + CC[0,2]adc 2 +CS[2,0]a 2 (i 
+DS[l > l]adt 2 + CS[0,2]ci 3 } 


(A5) 
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Appendix B 
Fourier Integration 
of 

Phase Function 


By defining that \j/j(t ) is the Fourier integration of phase function 


then 




2 k i(jk) 

f“ [1 

271 i(jk) 

— - — r [i 

2 7i i(jk) 


Pyfik) 2 + P^(ik) 


P 3j(ik) 2 + ik + P 4j 


] e ijkt ' d(jk) 


(Bl) 


i(jk) a ls 


i(jk) a 2j 


i(jk) + ja 3J i(jk) + ja 4 j 


] e’ jkt ' d(jk) 


By the characteristic of Fourier integral, the above equation can be inverted 
as 


[1 - PD,.] 
i(jk) 


-■if" V,(tO e- ijkt ' dt' 
2tc 

= i f" V*(tO e' ijkt ' dt' 
2tz 


(B2) 


ca 



The low limit of the integral has been changed from negative infinity to 
zero due to the fact that there is no negative time. By introducing the new 
dummy variable r such that 


r = i(jk) 


then eq. (B2) becomes 


[1 - PD/r)] 
r 


— f w,(tO e' rt ' dt' 
2n^ J 




(B3) 


Therefore, xj/^t’) is the inverse of the Laplace transform, i.e. 
Vj(t0 = SE- 1 ( [1 " PDj3 -} 


= S£- x {[- - 


±L_B 


r r + ja 3j r + ja 4j 
oV' 


= 1 - a y e' jv ' - a a e 


(B4) 
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Appendix C 
Newton’s Method 
for 

Finding Phase and Reduced Frequency 


As described in section 2.6, the governing equations for arbitrary 
motion are written as 


Fj = a - Oo cos(<|) + kt’) = 0 


F 2 = a + k sin(<j) + kt’) = 0 


Let the Jacobian matrix be defined as 



9F\ 

3k 

3<J) 

3F 2 

9F 2 

9k 

_ 


where 


(Cl) 


C.l 



-L- 7 = Oo sin(<{> + ktO, -^-1 
ok 9(j> 


=a 0 sin(4> + ktO 


3P 

— -i=a 0 sin(0 + ktO + a 0 k cos(<)) + ktO, — i=a 0 k cos(<{> + ktO 

^ d<|> 

The absolute value of the Jacobian becomes 


(C2) 


I Jacobian I 


- aFi aFa 
3k 9<J> 9k 


= -aosin 2 (0 + ktO 


(C3) 


By Newton’s method, the new values of unknowns k and 4> are 
calculated as 


Ki 


V 


'Ak; 


— 


+ 






A4> i_ 


(C4) 


where 


Akj 

A^ 


b - CF, f 2 ] 


Jacobian' 1 
I Jacobian I 


F,.^k - F,.^! 

9(f) 

-F + P ,.£l 

3k 2 3k 


I Jacobian I 


(C5) 
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The iteration will continue until both F x and F 2 reach 10‘ 6 
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